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Abstract 

When a thin elastic sheet crumples, the elastic energy condenses into a network of folding lines and 
point vertices. These folds and vertices have elastic energy densities much greater than the surrounding 
areas, and most of the work required to crumple the sheet is consumed in breaking the folding lines or 
"ridges". To understand crumpling it is then necessary to understand the strength of ridges. In this 
work, we consider the buckling of a single ridge under the action of inward forcing applied at its ends. We 
demonstrate a simple scaling relation for the response of the ridge to the force prior to buckling. We also 
show that the buckling instability depends only on the ratio of strain along the ridge to curvature across 
it. Numerically, we find for a wide range of boundary conditions that ridges buckle when our forcing 
has increased their elastic energy by 20% over their resting state value. We also observe a correlation 
between neighbor interactions and the location of initial buckling. Analytic arguments and numerical 
simulations are employed to prove these results. Implications for the strength of ridges as structural 
elements are discussed. 



1 Introduction 

The crumpling of a thin sheet is a phenomenon which we encounter every day, yet the equations governing 
crumpled systems are almost completely intractable without the introduction of drastic simplifying assump- 
tionsQ 0]. At the same time, this mundane occurrence exhibits some of the more intriguing behaviors of 
modern soft matter physics, such as phase transitions Q, scaling and energy condensation Q. 

For a large class of compressive boundary conditions, the energetically preferred configurations of crum- 
pled thin sheets consist of mostly flat regions bounded by straight folds and point-like vertices. Figure |l| 
shows an example of the resulting network of folds and points in a crumpled sheet of paper. One approach 
to analyzing such configurations is to treat them as patches of unstrained surface bounded by regions of 
discontinuous curvature |5j. Boundary layer solutions are then grafted to the regions of sharp curvature, and 
the total energies of the configurations are compared to find local or global energy minima. 

In the last several years the structure and energy of boundary layer solutions around straight folds and 
isolated vertices in crumpled sheets have been studied in detail, both from a physical perspective |^, 3|, [i . 



|, ||T^, |Tl|, m, |T||l|, |T|, |T|, 0, [T|, 1^, I2I, a nd mathematical perspective |, || pi ||, |25|, |26, 



27| . Related geometries such as thin film blistering |28, |2^, ^ |3^, thin viscous sheets |33| , thin film 



actuators |3J| , molecular sheets ^ , and the generalization of crumpling to higher dimensions |4|, |3^, 
|39t have also received attention. In particular, the boundary layer around a fold was extensively studied by 
Lobkovsky and Witten ||, 0, ^. They called the energetically preferred configuration a "stretching ridge," 
since it comes about through the balance of bending and stretching energy on the fold line, where both 
energies are of comparable magnitude. 

By viewing a crumpled sheet as a collection of ridges and vertices and adding up the known energetic 
cost of each unit we may arrive at a reasonable estimate of the total elastic energy in a crumpled sheet. 
However, part of the picture is still missing. We know from common experience that the crumpling of a piece 
of paper between the hands is a dynamic process, with the details of the final shape of the paper depending 
strongly on the history of applied forces and the effects of geometric frustration. In order to understand 

'^As Landau and Lifshitz dryly put it, (translated) |^ "These equations are very complicated, and cannot be solved exactly, 
even in very simple cases." 
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Figure 1: A typical crumpled sheet 
Image (a) shows a sheet of paper which has been lightly crumpled between the hands. Image (b) is the same 
sheet unfolded - lines and points resulting from plastic deformation show the former locations of folds and 
vertices in the crumpled state. Image courtesy the authors of p|. 



highly crumpled objects, which are clearly not free to find a global minimum of elastic energy, we need to 
know more about the energetic paths the membrane may take from one crumpled state to another. 

This thesis investigates the energetic pathway whereby one ridge buckles into several under the action of 
a compressive load. The work builds on Lobkovsky's scaling analysis of stretching ridges, though we differ in 
our treatment of applied forces at the tips of the ridge. Using improved simulational techniques and greater 
computing power, we investigate the buckling transition in far more detail than was previously possible. We 
then analyze the transition in a framework of stability and bifurcation theory, comparing and contrasting 
the transition on the boundary layer to the well studied subject of thin cylinder stability. 

We begin in Section ^ by reviewing the elastic theory of thin sheets, giving a brief derivation of the von 
Karman equations upon which our analysis is based. We then present Lobkovsky's derivation of ridge scaling, 
and his treatment of small perturbations to resting ridges. He chose a perturbative approach which assumed 
linear response to applied forces. After describing his method, we present an alternate approach which 
integrates applied forces into the original scaling equations. Our new approach provides better descriptive 
power when considering some special cases of applied forcing, since it does not assume that the applied forces 
are small. We then argue that our technique applies well to the highly relevant case of a ridge with inward 
point forces applied at its tips. Force scaling exponents are derived for this case. 

In Section ^ we present numerical evidence to support our scaling arguments. We devised a novel finite 
element program to increase the accuracy and efficiency with which we can simulate elastic sheets. This 
data was also presented in a companion paper fl^ . For most simulations presented here, the shape and 
boundary conditions of the grid were chosen to simulate a section of a cubical box, as shown in Figure |^. 
The grid covers one edge of the cube and has refiective boundary conditions. Simulating only this edge is 
equivalent to simulating a cube which is constrained to have 16- fold symmetry (all fold-lines equivalent). We 
argue that the reflective boundary conditions are representative of the real boundary conditions for a single 
ridge in a general crumpled sheet, since the forces which maintain the angle of a given ridge are most often 
exerted by the surrounding ridges. In later simulations we change the resting angle of the ridges, so that 
they do not correspond to edges of any closed polyhedral surface - the initial choice of the ^ dihedral angle 
corresponding to a cubic surface was arbitrary. Using data generated by these simulations, we demonstrate 
a scaling solution for ridges with inward forces applied at their ends. We also provide numerical evidence 
that the critical strain and curvature on the ridge at the buckling threshold scale with the same exponents 
as for ridges at rest. 

Finally, in Section ^ we consider the buckling transition. When the ridge is subject to strong enough 
forcing at its endpoints, it breaks into several ridges as shown in Figure We show that this transition is 
identical to the bifurcation by which thin cylinders buckle. We begin by reviewing the buckling transition of 
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Figure 2: Typical elastic sheet used in this study 
(a) The resting configuration of the simulated sheet with no external forces. It also shows the reflection 
planes to which the sheet edges are constrained, (b) How the simulated sheet is equivalent to one edge of 
a cube, when the mirror images of the sheet across the reflective planes are drawn in. The thickness of the 
sheet is .0004 of the the edge length and the Poisson ratio is 1/3. Darker shading represents higher strain 
energy density. The entire simulated sheet is uniformly darkened to distinguish it from its mirror images 
in (b). Slight numerical symmetry breaking between the left and right sides of the diamond created slight 
mismatches of the inferred surfaces on other faces, such as the right-hand face. The numerical grid is visible 
as a quilt-like texture. It has a finer scale at the edges and corners where curvature is larger. 



a cylinder under uniform axial compression. Thin cylinders subject to such forcing are observed to buckle 
in a regular diamond shaped pattern with azimuthal periodicity determined by their thickness and radius. 
The first bifurcation is shown to occur at a critical stress which is inversely proportional to the radius of the 
cylinder. We apply the relations derived for the cylinder to the geometry of the ridge boundary layer and 
show that it is consistent with the observed scaling of the ridge buckling transition. We also show that the 
small longitudinal curvature along the ridge and the non-uniformity of the curvature and strain on the ridge 
have only a weak effect on the buckling transition derived for the uniform cylinder. 

Adopting the hypothesis that the cylinder buckling mode accurately describes the buckling of stretching 
ridges, we make two previously untested predictions for the buckling behavior of ridges. Our first prediction 
concerns the normal mode which is associated with the buckling motion. We reanalyze our existing data 
from Section ^ to isolate the soft normal mode which becomes unstable as its associated spring constant 
passes through zero. Comparing the numerical results to the theory developed for the cylinder, we show that 
the rate at which the spring constant approaches zero with increasing stress along the ridge is on the same 
order of magnitude as for the normal mode associated with cylinder buckling. Our other prediction is that, 
for a given thickness to ridge length ratio, the buckling transition will occur at the same maximum ratio 
of stress to curvature. To show this, we simulate several variations on the ridge geometry described above. 
Although the different geometries buckle at different values of longitudinal stress and transverse curvature, 
they all buckle near the same value of the ratio of these parameters. 

We conclude in Section |^ by discussing the implications of this research for the strength of ridges as load 
bearing objects. The immediate consequence of our research is that ridges are not as strong as once thought. 
Since the ridge buckles when the total longitudinal stress along the ridge-line reaches a critical value, the 
pre-existing stress found in resting ridges makes them easier to break than stress free shells and cylinders 
with the same radius. Strategies for strengthening ridges are discussed, along with topics for future research 
concerning the buckling transition. 



2 Pre-buckling behavior of ridges 
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2.1 The von Karman equations 

We consider an idealized thin elastic sheet with uniform Young modulus Y. The sheet has constant thickness 
h which is much smaller than its spatial extent in the other two material directions. In the regime where 
elastic distortions are small and slowly varying on the scale of h, stresses along the thin direction can be 
neglected in comparison to those in the long directions. In this regime the thin direction can then be 
integrated out of the governing elastic equations so that the sheet is completely characterized by its 
2-dimensional center surface. 

We assume that our sheet has no intrinsic strain or curvature, so we may define on it material coordinates 
X C TZ^ ■ The embedding of the sheet into TZ'^ can then be expressed as some function r (x) of the material 
coordinates x. The strain is defined as the change in length element dl under the embedding: 



dl'^ = dl^ + 2jijdxidxj , (1) 



so that 



dr dr 

The most general quadratic form for the stretching elastic energy density in terms of the strain can be 
written as 

= 2(i^^t,2) i^^olv + Ve^k^Jl'y^J-fkl) , (3) 

where v is the Poisson ratio and eij is the antisymmetric tensor. 

The stress is defined as the variation of Cs with strain, dij = dCs/djij, so we can write £s — 
with 

r 

cr,j = ^ _ ^2 [lij + V e,k£jnki\ ■ (4) 

Because of the nonzero thickness of the sheet, there is also an energy cost associated with bending of the 
sheet out of its local material plane. Bending is quantified by the extrinsic curvature, which can be expressed 
as the component of the second derivative of f normal to the local material frame. 



dr 
dxidxj 



The inverses of the eigenvalues of the curvature tensor are the local principal radii of curvature of the sheet. 
Geometric constraints [^0[ Q require that the curvature tensor satisfy 

diCjk = dkCji (6) 

for a sheet with flat metric. 

The most general quadratic form for the energy density associated with bending can be written 

>Cb = - {nCijKCij + KGeikEjiCijCki) , (7) 

where the coefficient of kq is the familiar Gaussian curvature. 

To establish the connection between the curvature energy and the bulk elastic moduli, we introduce the 
so-called Monge coordinates. These coordinates locally parameterize the center surface of the sheet by 

f{x) = {x + u,y + v,w), (8) 

where x and y are the material coordinates and u, v and w are small deviations from the flat, unstrained 
state. 
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Microscopic considerations of the finite sheet thickness j|] yield an energy functional for 



r — 



24(1-1^2) 



(d.djw) + ve^kSji {d.djw) {dkdiw) 



(9) 



To lowest order in w, Eq.( |^) gives Cy — didjW, so we can immediately make the identification Cb = 
with 



KG 



12(1 -;/2)' 

Yh^v 
' 12(1-1/2)- 



(10) 



The equations of equilibrium can be found by setting the variation of the total energy of the sheet to 
zero. In the presence of an external pressure field P, the condition for an energy extremum becomes 



PSw = 



(11) 



Grouping terms for in and out-of-plane displacements yields the equilibrium conditions: 



djaij = 0. 



(12) 



Yh^ 



12(1 - Z.2) " ('^»^^^^) = ^' (13) 

Therefore the equilibrium condition can be written: 

nV^Cu = cFjkCjk + P- (14) 



Equations (^^ and (|1J) are not enough to completely specify the system, so another equation of state is 
necessary. An appropriate equation is Gauss' fundamental theorem of surfaces |4l| ], 

det Cy = didj-fij - V2tr7y , (15) 

which relates strain to Gaussian curvature. Together the force Eqs. |l^ and |l^ along with the constraint 
Eqs. ^ and ^ are enough to completely determine the equilibrium configurations of a sheet up to arbitrary 
translations and rotations. Equations ^ and |l5| are respectively called the force and geometric von Karman 
equations Q. 

As a consequence of Eq. ^ the curvature tensor can be written as the derivative of a continuous curvature 
potential 

a, = d.djf. (16) 

Here the potential f{x) is not identical to the local function w used above, but is approximately equal to it 
for nearly fiat surfaces. Also, Eq. ^ is automatically satisfied if we write <7ij in terms of a stress potential x, 

o'y = SikSjidkdix- (17) 

In terms of the potentials x ^^nd /, the von Karman equations assume the very compact form 

KV^f^[xJ] + P, (18) 

^V^X = 4[/,/], (19) 

where the brackets product represents 

[a, b] = {d^dpa) [d^d^h) , (20) 
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Figure 3: The stretching ridge boundary layer 
Images (a) and (b) show the geometry of a sheet with one sharp fold, (a) is a side view, while (b) is a cross 
section of the fold at the midpoint of the sheet. Images (c) and (d) show a representation of the "stretching 
ridge" configuration, in which the boundaries of the sheet are still required to make a sharp angle. Image 
(d) shows the same cross section which is shown in (b) - this image illustrates how the curvature across the 
fold-line is lessened in the ridge configuration, and the region of large curvature has a width on the order of 
the maximum radius of curvature. Image (c) shows how the softening of the curvature requires stretching 
of the sheet along its mid-line. The geometry shown here, used for the arguments presented in j^, requires 
extensional strain along the mid-line - other geometries can have a net compressive strain along the mid- line, 
but the driving balance between curvature and strain remains the same |7|. 
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2.2 Ridge scaling solution 

We wish to use the von Karman equations to study the the boundary layer around a folding line in an 
elastic sheet. This fold plus boundary layer configuration is what Lobkovsky et al. termed the "stretching 
ridge." j|, An intuitive picture of the structure of the boundary layer is presented in Figure |[ As we showed 
in the last section, for thin sheets the bending elastic modulus is less than the stretching modulus by a factor 
of h'^. Therefore, for very thin sheets with relatively free boundary conditions^ minimal energy configurations 
are mostly strain free, with large deformations concentrated around folding lines. Very near a folding line, 
the curvature approaches a scale where bending and stretching energies are once again comparable and the 
local configuration is determined by a balance between these two energies. 

Conceptually, as Figure]^ presents, the stretching ridge can be approached from the limit of zero thickness, 
where the curvature at the folding line becomes singular. The boundary conditions which create the fold are 
point-like vertices at points A and B, which are maintained at a sharp curvature. As the thickness of the 
sheet increases, it is energetically favorable for the middle section of the ridge to have a lower curvature at 
the expense of stretching along the length of the ridge. The boundary layer around the fold thus acquires a 
saddle- like shape as shown in Figure ||(c) and (d) . The width of the boundary layer is on the same order as 
the transverse radius of curvature R and is much less than the length of the ridge X. 

Now we proceed to apply the von Karman equations to the stretching ridge. This system has two well 
defined typical length scales ~ the sheet thickness h and the ridge length X. We can therefore rescale the 
von Karman equations into a more convenient dimensionless form by expressing all lengths in units of X 
and all energies in units of k. The von Karman equations then become 

A'v4x = -^[/"/]. (21) 

Here /, and P represent the stress potential, curvature potential and external normal forces in respective 
natural units of kX~^, X , and kX~^ . All derivatives are taken with respect to the dimensionless variables 
x/X and y/X. The dimensionless small parameter A is given by 



vWYTi _ 1 ( h 



X ^12(1-1/2) 



(22) 



We consider a sheet with edge boundary conditions sufficient to create a single ridge and no normal 
forces {P — Q). We define our coordinate frame so that the origin is at the center point of the ridge and the 
center-line of the ridge is parallel to the x material direction. Since A comes into the von Karman equations 
multiplying the stress source term, the possible configurations of a thin clastic sheet are well described by 
a stress free, A = 0, folding solution plus boundary layers at the fold lines. Lobkovsky's insight in was 
to try a scaling solution for the boundary layer of a single ridge which matched the / scaling of the outer, 
sharp fold solution. For a fold of dihedral angle tt — 2a across the line y = 0, f = a \y/X\. Accordingly, on 
the boundary layer / should scale with the same power of A as the dimensionless transverse coordinate y/X. 
He therefore tried a scaling solution of the form: 

/ = A'^/, X = A'x, y = x^y/x, i = x/x, (23) 

where the tildes denote dimensionless, scale free coordinates and functions. Taking /3 < gives the proper 
limiting case of sharp curvature at zero thickness. Substitution into the rescaled von Karman equations with 
P = and retention of only the largest terms yields 

A^'^l^-A'^-K,/], 

^''''-''^ = -l>^"[lf]- (24) 

^For a discussion of the role of boundary conditions, see [ ^ . 
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In the A ^ limit the two sides of each equation must scale identically. Solving for the exponents yields 

(25) 

This translates to X^/'^ scaling of the boundary layer width, A~^/'^ scaling of the transverse ridge curvature, 
and A^/"^ scaling of the strain along the ridge length. So, to within factors of order unity, the radius of 
curvature across the ridge is i? « XA^^'^, and the total width of the boundary layer should be about the 
same. Also, if we assume most of the elastic energy is concentrated on the boundary layer, a region of area 
X^/^X'^, then integration of the energy functionals in Eq. || and Eq. yields total bending and stretching 
energies which scale as k\~^/^ . 



2.3 Response to external forces 

In deriving the ridge scaling exponents, Lobkovsky et al. posited a set of "minimal" boundary conditions to 
create a ridge - purely normal forces are applied as necessary at the edges of the sheet to maintain straight, 
sharp folds. Relying on these boundary conditions as necessary, he derived functional forms for many of 
the geometric and energetic quantities of interest on the ridge |n, o]. However, the derivation of the generic 



scaling exponents presented in Section 2.2 does not rely on any specific boundary conditions, and so these 
exponents should apply to a much broader class of ridge configurations than just the "minimal" ridge. 

These scaling arguments still place a few limitations on generic ridges. An important ingredient in the 
ridge scaling derivation is that there are only two length scales, the thickness h and the ridge length X. 
This is the essence of our notion of an unperturbed, resting ridge - the boundary conditions do not impose 
another length scale on the problem. 

To quantify the response of ridges to external forces, we must study how the unperturbed ridge evolves 
as a new length scale is added to the system. Our common experience gazing at a crumpled sheet of paper 
tells us that real elastic ridges do not live in isolation, but are influenced by other ridges around them as 
well as by the global geometry of the sheet. Since ridge scaling is witnessed under such circumstances, the 
governing equations must be fairly insensitive to most small perturbations. 

Previous efforts Q to study the response of ridges to external forcing relied on a perturbative scheme 
which assumed linear response to small forces. The response of the ridge to in-plane external forcing at 
the edges of the sheet was deduced by first solving for the stress field such a force would produce in a 
flat sheet, then using this solution to modify the equations governing the ridge. A ridge scaling solution 
was then substituted into the modified von Karman equations and series expanded in the magnitude of the 
external force. This technique had the weakness that the resulting expressions for the force response had an 
undetermined exponent (though the possible values of the exponent were limited). The technique was also 
very complicated. 

We use a less specific, more intuitive approach to study the overall response of a ridge to external forcing. 
Instead of solving for the detailed behavior of the ridge strain and curvature as the force is applied, we focus 
on how the forcing must be rescaled to have an equivalent effect on ridges of different length ratios A. For 
example, forces applied normal to the sheet enter the von Karman equations through the term P in the 
force equation. If we add the dimensionless form of this term back to the ridge scaling form of the equation 
(Eq. we find 



A-^|^=A-Mx,/]+P, (26) 



where we have used the ridge scaling values for (5 and 5. In order for the solutions to / and x to remain 
scale invariant, the field P must obey the scaling form 

P = A-ip(5,y), (27) 

where P is a dimensionless. scale free function of x and y. Conversely, the existence of the scale-free function 
P gives a relation between two equivalent external forces on ridges of different A: 

(28) 
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Thus we assume that with proper rescahng of our applied forces, the fields / and % obey the same A 
scaling laws on forced ridges as on resting ridge. This notion allows us to make strong statements about 
the complex evolution of ridge shape with applied forcing. Still, the flexibility of our approach is greatly 
limited by the requirement implicit in Eq. ^ that the spatial dependence of the applied force also scales^ 
with A. A perturbation scheme with which we are not free to fix the location of our perturbing force seems 
to be of limited physical interest. However, this scheme is sufficient to study some special cases which are 
particularly important. For example, it is well suited to the problem of a ridge with external point forces 
applied at its vertices, since the spatial location of the equivalent forcing will clearly remain fixed as A is 
varied. The other benefit of this scheme is that it does not rely on assumptions of linear response, so it 
describes the force response over changes of order unity in the fields / and x- 



2.3.1 Point forces applied at ridge vertices 

External forcing applied to the sheet enters the von Karman equations via the term P or via boundary 
conditions at the sheet's edges. We described the proper rescaling of P in the last section - in this section 
we calculate the rescaling of a particular kind of in-plane forcing at the sheet boundaries. Here and in the 
remainder of this work we consider an external potential which applies point forces to both vertices at the 
ends of a ridge. Since the applied forces have only delta function spatial extents and are applied at points 
which remain stationary under ridge scaling, we do not expect them to destroy the spatial scaling of the ridge 
solution. Therefore we may reasonably expect to find that the equilibrium configuration of a ridge under a 
given compressive force is identical to rescaled configurations of ridges with different material thickness and 
properly rescaled external force magnitudes. 

To calculate the proper rescaling of the forces on the vertices for a similarity solution, we consider our 
forcing as a boundary condition consisting mainly of an in-plane point force. This force amounts to a point 
stress at the edge of the sheet with the form 

= Fo5{y)- (29) 

So, to find similar scaled configurations of the sheet, we must scale a^xx the same way Gxx scales on the ridge. 
Since ^xx scales as A^/"^, 

<Jxx Yh-ixx n\-^/^X-^xx, (30) 

where ^xx is a dimensionless, scale free function. To express a'^x in similar fashion, we first substitute the 
scale free y-variable y = A~^/^j//X, so that S{y) = X~^^^X~^S{y). Thus, the proper scale free force can be 
written in terms of Fo as 

Fo = ^Fo. (31) 

For reasons which will become clear in Section ||, we cannot measure the force applied to our ridge with 
very good accuracy. However, we can measure the inward displacement A of the ridge ends caused by this 
forcing. These two quantities may be related by assuming that the work done by equivalent rescaled forces, 
given approximately by -FqA, scales the same as the total energy of the resting ridge configuration. The 
total energy of a resting ridge scales as kX~^^^, so equivalent values of FqA/k will scale as A"^/"^. Given the 
scaling of Fo from Eq.^ the scale free A must be related to the actual displacement by 

A = A^/^XA. (32) 

This scaling result has the predictable implication that the macroscopic strain A/X scales identically to jxx, 
the local longitudinal strain on the ridge. 



3 Simulations 
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Figure 4: Finite elements for stretching and bending 
The strain on each triangle is computed from the change in relative positions of its vertices. Curvature 
on one triangle is computed from the relative heights, normal to the triangle surface, of the vertices of the 
triangle plus the additional vertices of its nearest neighbor triangles. 

3.1 Numerics 

We simulate an elastic sheet using a fixed triangular grid with variable grid spacing (see Figure |^). Strains 
and curvature are taken to be constant across the face of each triangle. Strain is calculated on each triangle 
by measuring the relative deviation of its vertices from their predefined strain free positions. Curvature is 
calculated on each triangle from the relative heights normal to the triangle surface of the three triangles that 
share sides v/ith it (see Figure The relative heights Zi of the six points (a) - (f) are fit to a function of 
the form 



ai + a2Xi + flaj/i + 04 {xif' + a^^x^yi + ag {yif' 



1,6 



(33) 



where Xi and yi are the material coordinates of the vertices, 
identification 



Curvatures follow immediately from the 



Cxx = 2 X 04, Cxy = (15, Cyy = 2 X ag. (34) 

These six simultaneous equations for the Zi in terms of aj are inverted at program initialization to save later 
computational time. 

The edges of the grid are constrained to lie in reflection planes of the minimal unit pictured in Figure |^. 
For calculations of the curvature on triangles bordering these planes, the triangles see mirror images of 
themselves across the planes. The mutual attitudes of the reflection planes are such that the sheet would 
meet them all at normal angles if it were perfectly flat except for one sharp 90° fold between two opposite 
corners. As shown in Figure ||, the single ridge with this geometry is conceptually identical to one edge of a 
cube surface, provided the surface is constrained to be symmetric. Two corners of the sheet are the vertices 
of the ridge that forms along the fold-line when the elastic energy is minimized. 

The gridding was determined in two steps: flrst the desired area was gridded with an equilateral triangular 

^This limitation is similar to Lobkovsky et al.'s requirement that his forcing be "compatible" with the boundary conditions 
for his "minimal" ridge. However, he still had complete freedom to study any form of boundary forcing in the material plane, 
since the minimal ridge boundary conditions only specify the normal forces. 
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Figure 5: Simulational geometry 
An equilibrium configuration of the simulated sheet (white grid) is shown for a thickness aspect ratio A of 
10~^. The daxk planes mark the location of the reflective planes to which the sheets edges are confined. 
Points A and B are the center-points of the repulsive potential used to press on the vertices of the 
simulated ridge. 
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lattice, then the positions of the grid points wore remapped by the simuhaneous transformation 



x' ^ f{x,y), y' ^ g{x,y), (35) 

where / and g are fifth order polynomials which are constrained to be stationary on the edges and along 
the midline of the grid. This mapping provided a fourth-order smooth gradient in grid spacings on the flat 
faces of the ridge while concentrating the lattice spacing at the vertices by a factor of 10"^ and across the 
ridge-line by a factor of 10^ compared to the flat regions far from the ridge. The concentrations factors were 
chosen arbitrarily, within the limits of the mapping, to make the gridding near the vertices as flne as possible, 
since this is the region of largest gradients in curvatures and strain. The gridding is visible in Figure ^. 

Bending and stretching energies were assigned to the curvature and strains on each triangle using the 
forms for elastic energy presented in |^ 

e_B = -^i^A{CiiCjj — eik€jiCijCki) (36) 
3 

es = —GA{ji^jjj +2j,jj^j) (37) 

where is the antisymmetric tensor, k is a bending modulus, G is the 2-dimensional Young's modulus, 
and A is the area of the triangle. The stress energy expression is appropriate for a material with Poisson 
ratio of 1/3. The coefficient of the Gaussian curvature energy is not consistent with that given by Eq. |l^ 
for a uniform elastic material, but was chosen to maintain consistency with Lobkovsky et al.'s simulations 
in |3| H The direct contribution of the Gaussian curvature to the bending energy is much less than 
that of the mean curvature. Separate simulations verified that changing the value of the Gaussian curvature 
coefficient had no discernable effect in our data. The physical thickness h of the sheet is equal to k/G. 

Pushing on the tips of the ridge is accomplished by introducing repulsive potentials of the form V{r) ~ 
Gp 1 1?'— Xn\ centered around two points, one on each line to which a vertex is constrained (points A and B 
in Figure |5|). The center points are located at a distance X from one another and symmetrically placed with 
respect to the middle of the ridge. These points lie where the vertices would be if the sheet were sharply 
creased - relaxation of the ridge curvature draws the vertices inward from these points for an unforced resting 
ridge. The benefit of this potential is that it acts mainly on a small but finite area around the vertex. In 
earlier simulations, simple pushing of the vertex itself led to local collapse of the vertex tip without imparting 
any force on the main part of the ridge. Cp was varied to apply different loading. 

An inverse gradient routine was used to minimize the total elastic and potential energy of the sheet 
as a function of the coordinates of all the lattice points for given parameters k, Y and Gp. 

Using this routine we found minimum energy configurations for ridges of aspect ratio A ranging from 
1.25 X 10^^ to 1.77 X 10^^. The upper bound on A was determined by the range of validity of the ridge 
scaling solution - above this value the width of the ridge becomes comparable to that the sheet. At the 
other extreme, for A < 10~^ the radius of curvature at the ridge line becomes comparable to the spacing of 
our lattice and the simulation ceases to be accurate. 

3.2 Findings 

The plot in Figure || shows scaling of the total elastic energy in the ridge versus A for ridges at rest (Cp = 0) 
and at the buckling threshold. The data shown here is for ridges with dihedral angle ^ . Scaling of the total 
elastic energy for the resting configuration is consistent with a kX'^I^ dependence, in agreement with prior 
theory and simulation 0, ^ . Figure || also shows that the elastic energy measured at the buckling threshold 
exhibits exactly the same scaling as on resting ridges. This suggests that the ridge scaling developed for 
the resting ridge is still applicable to the ridge with forces applied at its endpoints. As we predicted in 



Section 2.3.1 , this particular form of forcing potential should not destroy the length scaling of the ridge. 
The inset in Figure ^ shows that that the energy correction at the buckling threshold is nearly a constant 
fraction of the total ridge elastic energy. 

Scaling of the force response is verified by the existence of a similarity solution for the ridge shape as 
a function of tip displacement A (see Figure 0). We considered scaling of the equilibrium value of Gyy 
along a line in the material coordinates which bisects the simulated ridge line. As a consequence of the 
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Figure 6: Energy of ridges at rest and at the buckling threshold 

Straight lines are least squares fits to a scaling form y = ax^ . In this plot A ranges from 1.25 x 10""^ to 
1.77 X 10^''. The plot shows the total elastic energy [Eb + Es) in the sheet after minimization. The scaling 
exponent fit for the resting ridge values (lower line) was —0.32, the fit at the buckling threshold was —0.31. 
The inset shows the difference between threshold energy and resting energy in units of the resting energy. 
This energy ratio is best fit by a scaling exponent of 0.05 ± 0.02 and is consistent with a constant ratio. 



for an unforced ridge plots of Cyy x A^/^ versus y x A 



1/3 



along 



scaling exponents presented in Eq. 
this line should be independent of A. Extending this result to forced ridges, we found numerically that the 
rescaled cross-ridge curvature profiles were also identical for forced ridges with the same equilibrium value 
of A x A°'^^. This A rescaling exponent is very close to the theoretical value of 2/3 derived above. Plot (a) 
in Figure shows values of Cyy, along a line in the material coordinates which bisects the simulated ridge 
line, for several different sheet thicknesses and two different values of rescaled ridge tip displacement A. For 
comparison, the unsealed Cyy versus y for a particular rescaled A is shown in Figure |^(b). 

With our generic treatment of the ridge force response we can rescale the observed configuration at one 
thickness to that for the simultaneously rescaled thicknesses and applied forces. It must be noted, however, 
that scaling of the ridge response to forcing does not imply identical scaling of the buckling threshold. 
Buckling of the ridge signals a bifurcation in the allowed equilibrium configurations at a critical applied 
load This is a completely separate topic, which we treat in Section^. There, we introduce a model 

for the buckling transition which reproduces the observed scaling of the critical energy without assuming it 
a priori. 



3.2.1 Dihedral angle scaling 

In addition to the simulational data presented above for ridges with dihedral angle ^ , we also simulated the 
force response up to the buckling point for ridges with different angles. In this section, we present data for 
ridges with thickness aspect ratios A from 1.25 x lO^'^ to 2.5 x 10^'* and with dihedral angles from ^ to 
In each of these simulations, the sheet had equal side lengths as before and was held by reflective boundary 
conditions similar to those described above to form a ridge between two corners. 

Lobkovsky showed that for ridges with dihedral angle tt — 2a the elastic energy scales as kq^/^ . The 
best scaling fit to the resting energy of our ridges had an exponent of 2.31 (see Figure ||). Interestingly, for 
all dihedral angles and thicknesses studied, we found that the total elastic energy at the buckling threshold 
was always approximately 20% greater than the resting ridge energy. The constancy of this ratio is discussed 
in Section 



4.4.1 



once we have derived the buckling criterion. 
We also found for all cases studied that the total energy in the ridge as a function of tip displacement A 
was well fit by the quadratic functional form 



E = Eo 



^r(A-A„f , 



(38) 



where Eo is the resting ridge energy and Ao is a numerically fit zero offset. Typical values of Ao were found 
to be within 10% of the resting ridge zero offset. The value of F was found to be nearly independent of 
dihedral angle, but was well fit by the scaling form T « 3.2 A^^-^^ k/X"^ (see Figure ^. For comparison. 
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Figure 7: Similarity solution for the ridge response to forcing 
Both plots show Cyy, the curvature across the ridge-line, versus the y material coordinate on the line which 
bisects the ridge-line. The data is for sheets with seven different values of A, ranging from 1.25 x 10~^ to 
1.25 X 10"''. Plot (a) shows Cyy x (X/Xi)^^'^ vs. y x (X/Xi)^^^'^ for the ridges at rest and for ridges with 
inward vertex displacement A(A) = A^*^' x (A/Ai)" '^'^, where A is measured from the resting vertex positions, 
Ai is the aspect ratio for the thickest sheet, and A^°^ was vertex displacement at the buckling threshold for 
the thickest sheet. The profiles with the large central peak are the buckling threshold values. (The small 
dimple in the data at y = is a numerical artifact due to a discontinuity in the gridding density across 
the ridge-line. For finer gridding this dimple goes away, while all other local values of curvature remain 
constant.) Plot (b) shows unsealed Cyy versus y for the buckling threshold profiles plotted in (a). 
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Figure 8: Scaling of ridge energy with a 

Data is shown for ridges with thickness aspect ratio A ranging from 1.25 x 10~^ to 2.5 x 10"''. The dihedral 

angles of the ridges are given by tt — 2a. The energy values for each thickness were rescaled by the predicted 
energy scaling factor A^''"^ so they all lie on a common line. Points on the lower line were for ridges at rest 
and points on the upper line were for ridges at the buckling threshold. The lines are fits to the data for 
thickness A = 5.6 x 10"'*. The a scaling exponents for this thickness were 2.31 for resting ridges and 2.29 
for ridges at their buckling threshold. The thinnest sheets ("V" symbols) do show some deviation from the 
top line fit. 
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Figure 9: Scaling of F with a 

The compressibility modulus T was calculated for ridges with thickness aspect ratios A from 1.25 x lO^"' to 
1.25 X 10"* and with dihedral angles from ^ to The line is a scaling fit for data with dihedral angle ^ 
and has an exponent of —1.65. 
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Figure 10: Diamond shaped buckling mode for a thin cyhnder 
This image was taken from . On the right is a photograph of an aluminum cylinder buckled by application 
of straight downward forces at its ends. On the left is a relief map of a simulated sheet with the same buckling 
pattern. 



the elastic energy of a thin strip of length X and width w whose ends are compressed inward by length A 
is approximately given by i? « ^Yh{/^/ X)'^Xw. If we take the width to be the ridge width w = X\^^^ and 
substitute Yh = k/K^ = k\~^ /X^, then the energy becomes E w iA~^/'^/c/X^A^. Thus the compressibility 
of a thin flat strip with width on the order of the ridge width is also of order A^^/'^k/X^. 

4 The buckling transition 

Detailed study of the buckling transition is complicated by an apparent discontinuity between pre-buckled 
and post-buckled states. When a sheet of paper or a tin can is buckled new ridges appear suddenly, often 
accompanied by the popping sound of energy release p9[ . In simulations we find that ridges buckle at 
a repeatable value of inward tip displacement, but immediately after the buckling transition the minimal 
energy configuration of the sheet contains a fully formed new ridge. This is counter to our intuition that the 
transition should be continuous Q, in which case the immediately post-buckled state would be only 
infinitesimally different than the pre-buckled. More importantly, our inability to observe intermediate stages 
in the growth of the buckled state prevents us from directly seeing the shape of the assumed normal mode 
against which the ridge becomes unstable. 

In this section we study the buckling transition in greater depth, through more detailed analysis of the 
simulational results presented earlier. The consequences of our analysis also lead us to run further, more 
specialized simulations. Our aim is to tie the observed behavior at the ridge buckling threshold to well known 
results concerning the buckling of thin elastic cylinders. 

We begin by reviewing the salient features of thin cylinder buckling. Under uniform axial compression 
applied at its ends, a thin elastic cylinder will pass from a state of uniform uniaxial curvature to an axially 
periodic diamond shaped buckling pattern (see Figure [lC|) . Although the details of the geometry are different, 
we show that the buckling of elastic cylinders or sections of a cylinder is determined by elastic terms analogous 
to those which dominate the behavior of stretching ridges - namely the cross-ridge curvature and the ridge- 
line strain. We proceed to apply the analysis developed for the cylinder to elastic ridges, finding that the 
scaling it predicts for the buckling transition is consistent with our simulational observations. We then do a 
normal mode analysis of the buckling transition to determine the scaling of the lowest mode as it approaches 
the point where the single ridge configuration becomes unstable. We numerically compute the lowest several 
normal modes for ridges under various degrees of compression, and show the appearance of a mode with 
swiftly decreasing eigenvalue which we believe accounts for the buckling instability. Quantitative aspects of 
this mode's shape and its approach to zero eigenvalue are discussed. 

Finally, we argue that the apparent differences between the elastic cylinder and the stretching ridge are 
inconsequential in regard to the buckling instability. The main distinguishing trait of the stretching ridge is 
that it maintains high internal stresses even when it is at rest. We show that the critical load at which the 
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Figure 11: Coordinates in cylinder frame 



ridge buckles is determined by the ridge curvature and the total strain along the length of the ridge, including 
the pre-existing stress and the additional stress resulting from the applied load. To demonstrate that the 
buckling transition only depends on two parameters, the transverse curvature and the total longitudinal 
stress, we present data from several additional simulations with different geometries. Though these new 
simulated sheets buckle at different total stresses and curvatures, they all buckle at the same ratio of these 
two quantities. We also discuss an observed universality of the additional energy required to break a ridge. 

As a last note, we address the jump in position and energy at the buckling transition. Even though the 
appearance of unstable modes for cylinders is the result of a continuous bifurcation in phase space, cylinders 
also jump discontinuously in energy upon buckling. We justify this jump in terms of the non-linear growth 
of the buckling mode. We also argue that the final wavelength of the buckling pattern on a ridge need not 
be the wavelength of the instability. 

4.1 Stability of thin elastic cylinders 

Our treatment of the stability of thin elastic cylinders and sections of a cylinder under a compressive load 
mainly follows that presented in |50[| . To avoid additional boundary conditions, we make our argument for 
a complete cylinder - however, the stability condition we find is local, so it can be directly applied to an 
angular section of a cylinder with the same local stress and curvature fields. We discuss local buckling and 
angular shell sections at the end of this section. 

We consider an elastic cylinder of thickness h, length L, radius R and Young's modulus Y. To assess the 
stability of the cylinder under a compressive force F = 27: Ra uniformly applied along its edges, we consider 
the stability of the force von Karman equation, Eq. ^ against infinitesimal displacements. The force F is 
conveniently expressed so that the resulting longitudinal stress in the unbuckled cylinder is a. Under the 
application of the force at its ends, the cylinder will naturally undergo some compression along its length, 
which will in turn cause it to expand radially. We are not concerned with these distortions, but only consider 
them as the equilibrium solution to the von Karman equations to which we add an infinitesimal displacement 
which will grow into a buckled solution. We define a local coordinate system everywhere on the cylinder 
with X direction along its length, y azimuthal, and z normal to the surface. Infinitesimal displacements in 
these three directions are labeled u, v, and w respectively. In this frame, z points radially inwards. 

Because our cylinder already has a curvature field Cyy — l/R, the relation between the additional 
displacements u, v, and w and the resulting strain and curvature fields is not as simple as for flat sheets. Up 
to an additive constant, the local embedding of the cylinder into TZ^ is given by: 



This expression accounts for the rotation of our local frame, as demonstrated in Figure |l^. Referring to Equa- 
tions I and |, the expression for the additional strain and curvature due to our infinitesimal displacements 
are, to first order in u, v and w, 




(39) 




(40) 
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^■^ dx^' dy^ Rdy' dxdy Rdx' ^ ' 

Here the primes denote the infinitesimal corrections to the equiUbrium fields. For equilibrium small dis- 
placements, the terms in Eq. ^ involving derivatives of v are typically much smaller than those involving 
derivatives of w, so they are neglected in the following treatment. The solutions we find for u, v and w are 
consistent with this approximation for the values of h/R considered. 

Empirically, thin elastic cylinders typically buckle in a diamond pattern like that shown in Figure |l^. 
The longitudinal and azimuthal periodicities of these patterns vary. To test the stability of the shell against 
these buckling modes, we consider a family of infinitesimal displacement of the form: 

w = yle^"/-Rcos(7V2//i?), 
V = Be^'^^/^smiNy/R), 

w = Ce*"/-«cos(iVj//i?), (42) 

where the periodicities r and N are free variables. This form of the buckling mode neglects boundary effects, 
and so it is most accurate for very long cylinder, for which R/r 3> L. 

We neglect higher order corrections due to the geometric von Karman equation in our consideration of 
infinitesimal displacements, but the force von Karman equation is not valid unless the additional constraint 
of in-plane force equilibrium is satisfied. From Eq. |l^, this requires dicrlj = 0. Using the relation of stress to 
strain from Eq. ^ combined with Eq. this yields two equations: 



l±^d,dyU + ^dlv + - ^dyW = 0. (43) 

These equations can be used to solve for A and B in terms of C, yielding: 

A vr^-N^ B ^i2 + u)r^ + N^ 

= —ir = N- - (44) 

C (r2-FiV2)2' C (r2-hA^2)2 ■ ^ > 

The equilibrium cylinder configuration becomes unstable to the combined displacements in Eq.^ when 
the resulting normal force, calculated from the force von Karman equation, goes to zero. To linear order in 
the small displacements, this condition will be met when the terms linear in A, B and C satisfy 

«:V2C^ = -la;, + (TC;,. (45) 

Using Eq. ^ to express a'yy in terms of and 7^^ , the above equation can be written in terms of the small 
displacements as: 

n^\dlw + dlw) + adlw - ^ (dyv + vd.u - ^) = 0, (46) 

where we have used Eq. ^ to express all elastic moduli in terms of k. 

Substituting the buckling form of Eq. ^ directly in Eq. ^ and using Eq. 44 to eliminate A, B, and C 
yields the bifurcation condition 



/l2 

If we define 

„2 



r 

(r2 + 7V2)2 ' 
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Figure 12: Relation between r, N, and 77 
The solid line is the higher root of Eq. pA, while the dashed line is the lower root. 



we can write Eq. ^ as 



Solving for a yields 



12(1 



-1]^ - aR^T] + K = 0. 



(l2(l-z.^)(f)%^ + l)^ 



This function has a minimum in 77 when 



77 = ,7,, = (12(1 - 7.2)) 



2.^-1/2 h 



R 



(49) 



(50) 



(51) 



Where 77c/ is called the "classical" buckling value of 77. The corresponding minimum value of a is denoted as 



(52) 



This is referred to as the "classical" value of the buckling stress The classical stress is the smallest 

applied load under which the cylinder can buckle, provided that a solution for the corresponding value of 7/ 
is allowed.^ 

Our limit of a thin sheet corresponds to the limit 77c/ — » 0. The variables r and N are not uniquely 
determined by this stability treatment, beyond the requirement that they satisfy Eq. ^ However, for each 
set of r and N there is a unique value of 77 and therefore of cr > ad , determined from Eq. at which the 
cylinder becomes unstable to buckling with that mode. If r and N are taken to be continuum variables, 
then there is an entire family of solutions which satisfy t] — rjd. In real cylinders the preferred buckling 
wavenumbers are determined by boundary conditions and initial imperfections. Azimuthal periodicity and 
the requirement that the wavelength be commensurate with the finite length of real cylinders severely limits 
the number of allowed solutions for cylinders which are relatively thick. In such cases, there may only be 
a handful of allowed combinations of r and N for which rj « 7/c/, with corresponding threshold instability 
value of a close to a^. Also, it has been shown |pT| that initial curvature imperfections in the cylinder can 
break the degeneracy in r and N . 



*It has been found experimentally that, depending on boundary conditions, 
diamond pattern near the cylinder ends can lower the bucking stress by up to 50% 



tial imperfections and deviations from the 
pl| ]. However, our treatment is sufficiently 
accurate to establish the length and energy scale at which cylinders will buckle when boundary effects are neglected. 
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Without yet placing any constraints on the allowed values of r and N ^ we explore the limits imposed by 
Eqs. and pll. Solving Eq. Bq for the wavcnumber r in terms of N and r] gives 



1 



r^=^\l±{l-mSy'^\. (53) 

Thus r can take values between zero and rj~^/'^ and N can range from zero to ^77"^^^. Possible combinations 
of r^/rj and N ,Jrj are shown in Figure |2|. For most allowed combinations of r and iV, r is on the order of 
7^-1/2^ For small r and N-q^/'^ we can expand the lower root of Eq. finding r w ■q^/'^N'^ . Substituting the 
value of rjci from Eq. |l] gives two limits for large R/h: 



^/R/h > r > y/h/R, (54) 

where for the lower limit N <^ ^ R/h. Implications of both bounds for ridge buckling will be described in 
later sections. 



For comparison with later results, we note that the classical breaking stress given in Eq. 52 implies a 
breaking strain of order 

where we have assumed ayy « and have used ^' = 1/3 to correspond with our simulations. Thus, the total 
energy input required to buckle the cylinder is 



If , , r,, kL 



(56) 



which is on the same order as the energy required to bend the cylinder out of a flat sheet 

\jnCldA = .'^. (57) 

As we mentioned above, the threshold condition for cylinder buckling can also be considered as a local 
condition on an angular section of a cylinder. The form of the buckling displacements which we posited 
in Eq. |4^ is a global motion, but we can construct a localized wave packet from combinations of these 
modes. Especially for short wavelength modes, there are many different modes with corresponding rj near 
7]ci, SO a localized packet can be constructed with critical stress very close to aa- The only requirement for 
this to work is that the width of the packet be at least a couple times longer than the principal buckling 
wavelength. From Eq. |5j, the smallest packet must therefore have a width of at least 2'KR/r = 2'Ky/hR. All 
of the buckling motion may be localized within such a self-contained wave packet, so the buckling threshold 
condition should be determined only by aspects of the cylinder defined within the wave packet. In this 
instance, it is appropriate to rewrite Eq. ^ with the local curvature Cyy in place of so the local 

buckling condition for v = 1/3 is: 

l.x/Cyy = 0.61/1. (58) 

Thus, if the ratio of strain to curvature locally surpasses the threshold value in Eq. ^ over a region of spatial 
extent greater than 2'Ky/hR, local buckling can occur just at that point. Experimentally, cylinders are often 
observed to begin buckling locally instead of all at once, due to inhomogeneities in the cylinder material and 
uneven forcing |5l|. 

4.2 Application to ridge stability 

It is immediately apparent that the elastic terms which dominate the buckling behavior of cylinders are anal- 
ogous to those which determine the scaling behavior of ridges. The largest terms in the stability condition. 



Eq. 45, are those proportional to the transverse curvature 1/i? and the longitudinal strain a. Likewise for 
the ridge, the balance between transverse curvature and longitudinal strain determine the ridge's shape and 
energetics both at rest and under compression. Furthermore, if we substitute the scaling form of the ridge 
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Figure 13: Distorted cylinder with one large and small curvature 



curvature l/i? = Cyy ~ A ^^"^ jX into Eq. p9 for the critical strain of the cylinder under applied load, we 



find ^ A^/-^, which is exactly the scaling we observed for the ridge in Section 3.2. 

We therefore anticipate that the buckling mode of the ridge should be roughly the same as that of a 
cylinder of material thickness /i, length L — and radius R = A^/"^X, where X is the ridge length and A is 



the thickness aspect ratio defined in Sec. 2.1. We picture the boundary layer of the ridge as behaving like 
an angular section of a cylinder, with some semi-rigid boundary conditions at the edges where the boundary 
layer meets the ridge flanks. This picture makes two assumptions: first that the real buckling mode of the 
ridge is localized on the boundary layer and second that the longitudinal curvature on the ridge-line doesn't 
change the critical strain scaling. The first assumption is supported by numerical evidence in the next 
section. Also, the ridge-line is observed to absorb nearly all the stress of our applied load without noticeably 
changing the shape of the ridge flanks, so it is plausible that the boundary layer only sees the rest of the 
sheet as a set of boundary conditions. 

In the remainder of this section, we calculate the correction to the buckling stress resulting from the 
non-zero longitudinal curvature. For ridges with typical aspect ratio A = 10^^ or less we anticipate that 
the longitudinal curvature is too small to have a pronounced effect on the buckling transition. We argued 
earlier that the longitudinal radius of curvature goes to zero as y^l'^X^ so for an aspect ratio of A = lO^'^ 
this radius of curvature already is of order lOX. The buckling wavelength couldn't be longer than the ridge, 
so the buckling deformational mode should not be strongly affected by the smaller curvature. 

In order to treat the longitudinal curvature rigorously, we repeat the derivation of the last section for 
a surface with curvature R\ in the y material direction and — -R2 hi the x direction. We take i?2 ^ R\- 
These curvature fields cannot globally describe a real surface, but for Ri small, we may picture a surface 
such as that in Figure where R\ is nearly constant over the length of the object. For such a surface, the 
strain-displacement relations analogous to Eq. ^ become 

, du w I dv w I \ f du dv\ 

^''■'^Ii^R'2^yy^d^~R'i^''y^2\d^^Ii) ^ ' 

and the in-plane strain equilibrium equations become 

dlu + ^—^dlu + ^—^dxdyv - -^dxw + -^d^w = 0, 
I A Hi H2 

^d^dyU + ^dlv + dlv - ^^dyW + ^dyW = 0. (60) 



If we substitute the buckling mode from Eq. 42 into the above equations, the ratios of buckling coefficients 
become 

A —ir ( 2 Ar2 ^1 



C (r2+iV2) 



C (.2+^2)2v(2 + -y+^^-^(-^^-^))- (61) 
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Eq. ^ for force balance of the infinitesimal displacement gains several additional terms from the new 
curvature, 



kV w + ad^w - — — dyV + vd^u - — + v— 
Kill \ Ki 1X2 



Solving as before yields the equation 



12k „ WW. , . 

^U. + .5,.-.- + -H0. (62) 



/l2 

If we define 



12(1 i^')i?t«r^_^\%(^2^^2)2(^(^2^^2)2_^^y)^0_ (g3) 

Hi K2 J 



^' = (^2+^2)2' ^' ^2_|i^2' ("^^^ 

We can proceed as before to find the lowest allowed value for a' as a function of ry'. The corresponding value 
for cr is 

J.2 _ Rl]\f2 

^ = ^VW^)^^^,. (65) 

This stress is ^ower than the breaking stress for a cylinder, but will approach the same value for |^ ^ 0. This 
justifies our above assumption that the minor Cxx curvature has a weak effect on the buckling threshold. 



4.3 Numerical investigation of buckling modes 

As mentioned above, our initial observations of the buckling transition came from detailed simulations in 
which we approached the buckling threshold with very small steps in inward vertex displacement. In all cases, 
the buckling transition was accompanied by a downward jump in total ridge energy and a large decrease 
in the strain along the ridge. There was no noticeable change in the vertex position after buckling. This 
assures us that the sudden jump is not aided by any work done on the ridge due to the springiness of our 
potential. 

The transition to the buckled state was of course accompanied by the appearance of additional ridges 
and vertices in our simulated sheets. Since the new ridges and vertices appeared in regions which were not 
finely gridded on our lattice they were often accompanied by curvatures which were large on the local scale 
of inverse lattice spacing. For this reason we do not claim that the buckling patterns we observed match 
in detail the real buckling pattern of a physical sheet, though they should be qualitatively correct. The 
types of buckling pattern we observed are illustrated in Figure ^ For relatively thick sheets, the initial 
buckled state was like that in Figure ^(a) - one new ridge appeared across the original ridge, with length 
on the order of the unbuckled ridge width, positioned about two fifths of the way along the original ridge. 
In simulations with the same sheet thickness but larger steps in hard-wall potential position (so that the 
initial step across the buckling threshold pushed deeper into the buckled state), the first observed buckled 
state was like that in Figure ^(b) , with two new ridges positioned symmetrically about the midpoint of the 
original ridge. For thinner sheets like that shown in Figure ^(c) the buckling pattern consisted of a number 
of smaller ridges. As with the larger ridges in (a) and (b), these smaller ridges had length on order of the 
original ridge width and were clustered near locations about two fifths of the way along the original ridge. 
This is consistent with the buckling patterns on a tetrahedron seen in In that system, which is closely 
related to our own, buckling was accompanied by the appearance of two large transverse ridges, which were 
symmetrically spaced about one quarter of the way along the original ridge. In only highly buckled 
states of the tetrahedron were observed, so we cannot be sure of where the new ridges initially formed on 
the tetrahedron (we have observed our o wn ri dges to change position slightly with the growth of the buckled 



state). This issue is revisited in Section 4.4.2 
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Figure 14: Observed post-buckling configurations 

Images (a) and (b) show the buckhng pattern for ridges with thickness aspect ratio A — 1.25 x 10^'^. The 
configuration shown in (a) was at the smallest step past the buckling threshold that we simulated for this 
thickness. Image (b) shows a ridge with tip displacement further past the threshold value. Image (c) shows 
a buckled configuration for a ridge with thickness aspect ratio A = 2 x 10~*. Lighting and shading were 
chosen to emphasize physical features. 
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4.3.1 Normal modes 



Since there seems to be some randomness involved in the selection of the postbuckled state accessible to 
our simulations, and since directly after buckling the system passes through non-equilibrium states of inter- 
mediate energy which we cannot directly observe, we desire to learn what we can about the buckling mode 
before the ridge buckles. Prior to buckling, the infinitesimal displacements given in Eq. ^ would result in 
a restoring force normal to the surface opposing the growth of buckling mode. As the stress a approaches 
the buckling stress CTcZ, this restoring force goes to zero. Because the boundary conditions of a cylinder or 
a ridge enforce selection rules on the allowed buckling wavelengths, the buckling mode must be part of a 
discrete spectrum of eigenmodes for motion of the sheet. Very near the transition, the buckling mode is very 
soft, and at some point before the associated eigenvalue goes to zero it must have the lowest normal mode 
eigenvalue. 



Ignoring the small corrections calculated in Section 4.2, we can surmise how quickly the buckling mode 
approaches zero eigenvalue with applied stress by expanding Eq. ^ around a — (ad — crs). For non-zero as, 
the left hand side of Eq. ^ will not equal zero, but will instead equal the restoring force per unit area linear 
in the buckling mode amplitude, which we denote by Pg. If we also ignore the small changes in transverse 
curvature 1/R with the applied stress, then the ad term cancels all the terms on the left hand side of Eq. ^ 
except the as term, and we are left with 

Ps = asdlw = —J—e^rx/R cos{Ny/R), (66) 

where we have substituted the form for the buckling mode from Eq. The work required to cause this 
displacement is then 

Ws = J Psw^\Psw\RX, (67) 

where R is the ridge radius of curvature, X is its length, and we again assume that the buckling motion is 
confined to the ridge boundary layer. The resulting expression is quadratic in the displacement amplitude 
C, and so by our previous ansatz we can identify it with the spring constant of the eigenmode which leads 
to buckling, by Ws = \KsC^ ■ Expressing as in units of ad, we can write 

as = as /ad (69) 

Finally, if we substitute in the scaling form of the ridge curvature for 1/R and the limiting scalings of the 
classical buckling value of wavevector r, the expected spring constant of the buckling mode becomes 



-^\-^/^ds r ^ ./R/h, 
j^N*\~^ds r - y^hjRN^, N < y^R/h. 



Ks^iTL.^^^i , „ ^ (70) 



Clearly, for A < lO^'^ the factor of A^^/"^ will become very large, so the shortest wavelength buckling modes 
will approach zero very quickly on the scale of the ridge parameters. Even the lower value of X~^N'^ should 
be a pronounced feature in the normal mode spectrum for any TV of order unity. 

We can also calculate the possible wavelengths of the buckling modes on the ridge- line by substituting the 
scaling dependencies into the two limits of r presented in Eq. |5^. The high wavenumber cylinder buckling 
mode has r ~ \/R/h ~ A~^/^, while the low wavenumber mode has r ~ y/h/RN'^ ~ A^^'^A^^. Thus, the 
limiting wavelengths of the classical buckling mode are 

C ^ 27ri? ^ h^AA^/a r ^ ^/W/h, 

r |27rAiV-2 r ./h/RN^ , N <^ y/W/h. ^ ' 

So for an aspect ratio of A = 10"'^ the shortest buckling wavelength is of order X/20. The longer wavelength 
is independent of A. On the ridge, the transverse wavenumber TV does not have a lower bound, but the 



24 



(a) 300 



X 250 - 




I 1 1 1 Li 1 

0.05 0.1 0.15 0.2 0.25 

inward tip displacement (X) 

(b) 90 |- , 1 , -. , 1 




I 1 1 1 1 : 1 

0.20645 0.20655 0.20665 0.20675 0.20685 0.20695 
inward tip displacement (X) 

Figure 15: Evolution of modes 
Both graphs plot values of the effective spring constants Kg as a function of ridge tip displacement for 
eigenmodes of a ridge with aspect ratio A = 2 x 10~^. The ridge tip displacement is nearly linear in the ridge 
strain and stress, so for these graphs ag ~ (0.2 — A)/0.2. The top graph plots the eight lowest eigenvalues 
at several different ridge tip displacements (under application of inward external forces at the tips). The 
dashed line in (a) has a slope of approximately 4 x 10^ x'^a-s ' lower graph is a closeup of the lowest four 
eigenvalues very near the buckling threshold. The dashed line in (b) has a slope of approximately 10^ x^g^ 

(for comparison to scaling values, (0.002)"^ = 5 x 10^ and (0.002)"'^''^ = 2 x 10**). The labeled points (A) 
and (B) correspond to the modes pictured in Figure |l6|(a) and (b) respectively. 



longitudinal wavelength (ci will presumably be a half-integer fraction of the ridge length. We therefore 
expect that the longest of these buckling wavelengths will be on the order of the ridge length. For — X , 
the corresponding value of is iV = V27r w 2.5. 

Numerically we looked for the buckling mode among the lowest normal modes of our simulated sheets 
for equilibrium configurations from zero forcing up to the buckling threshold. We found the modes by 
analytically calculating the the matrix of second derivatives of the total elastic energy for the equilibrium 
positions of the sheet, then using a block-Lanczos algorithm to find several of the lowest eigenmodes of 
this matrix. The eigenvalue corresponding to the buckling mode is precisely the spring constant Ks defined 
above. We used the "Underwood" implementation of the block-Lanczos algorithm, which is freely available 
on the NetLib online archive |5^. The block-Lanczos method is efficient at finding extremal eigenvalues 
and eigenvectors for large sparse matrices - our matrices were large by virtue of the large lattice size, but 
sparse since local curvature and strain fields at any lattice point are determined by relative positions of other 
grid points only up to a distance of next-next-nearest neighbors. The numerical values of elements in our 
second derivative matrices for very thin sheets differed by up to four orders of magnitude, so convergence 
of the Underwood routine was slow, taking up to several days to recover eight eigenmodes on a 700MHz 
Linux-based computer. 
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Figure 16: Representative eigenmodes 
In both images, the x and y coordinates are the material coordinates of the sheet, while the z coordinate 
is the eigenmode motion normal to the ridge surface. In these sheets the ridge-line extends from the upper 
left to lower right corners. The top image shows a longer wavelength mode which covers the length of the 
ridge. The the eigenvalue and tip displacement for this mode are label by point (A) in Figure ha(a). The 
bottom image shows a short wavelength mode. The the eigenvalue and tip displacement for this mode are 
label by point (B) in Figure ^5|(b) . 



26 



Figure |5| shows the evolution of the eight lowest eigenvalues as a function of vertex displacement for a 
ridge with aspect ratio A = 2 x 10"'^. The eigenvalue evolution is qualitatively the same for thinner sheets as 
well. Over a large range of inward vertex displacement, the lowest modes all have nearly constant eigenvalues. 
The modes contain an assortment of motions which are either global or localized on the boundary layer or 
the ridge flanks. As the vertex displacement (and therefore the ridge stress) is increased, the eigenvalues 
corresponding to long wavelength modes localized on the ridge begin to drop more steeply. For example, the 
strongly sloped line in the upper right corner of Figure p^(a) corresponds to the mode shown in Figure p^a). 
The eigenvalue for this mode has a slope on the order of —4 x 10^ (in units of ^/a^), which is just one order 
greater than the minimum slope predicted by Eqs. and ^ for the \~^N scaling associated with long 
modes of wavelength X {N = V^n). Similar modes, with wavelengths 2X/3, X/2, 2X/5, etc., were found 
higher in the eigenmode spectrum. Near the buckling threshold, these modes were also seen to approach 
zero eigenvalue with slopes greater than, but on the order of, that for the wavelength X slope described 
above. The computational time required to calculate higher modes prevented us from studying them in 
greater detail, but they behaved fundamentally the same as the mode shown in Figure |l^(a) . 

None of these modes ever reaches zero value before the ridge buckles, however - instead a very localized, 
short wavelength mode like that pictured in Figure |l^(b) appears suddenly, with very sharply dropping 
eigenvalue, just before the ridge buckles. The evolution of the eigenvalue associated with this mode is shown 
in Figure ^(b). The short wavelength mode has a wavelength on the order of the lattice spacing and is 
asymmetric about the center point of the ridge. For ridges with aspect ratio A = 2 x 10"'^, Eq. ^ gives 
a minimum wavelength of approximately X/IO, which is on the same order as the mid-ridge local lattice 
spacing. Thus the observed short wavelength mode is on the order of the minimum allowed wavelength, 
and should therefore have a spring constant near that predicted by Eq. ^ for X^'^^^ scaling of Kg. As 
Figure |l^(b) shows, the final slope of the eigenvalue as it approaches zero is indeed on the right order of 
magnitude to fit the cylinder buckling theory. 

Both the long and short wavelength modes seem to obey the scaling of cylinder buckling modes as they 
approach instability. That the short wavelength mode reaches zero eigenvalue first in every case could be 
due to some suppression of the long mode, either by boundary conditions or by the changing geometry of the 
boundary layer along its length. Also, the short wavelength mode may well be enhanced by lattice effects, 
and therefore should be more prone to cause buckling. The short wavelength mode is also enhanced by its 
high localization - since the stress and curvature are not uniform along the ridge (see Figure ^7|), localized 
patches of the ridge-line will meet the stress to curvature threshold criterion, Eq. M, before it is satisfied 
globally. In any case, the theory developed for cylinders at the beginning of Section^ allows for both these 
families of modes approaching zero spring constant just before the ridge buckles. The observations of these 
modes and their matching eigenvalue slopes strengthens the connection between ridge and cylinder buckling. 

The spatial extent along the ridge of the short wavelength mode envelope was observed to be independent 
of the sheet thickness. The longitudinal wavelength of the buckling was at the lattice spacing. Figure |l7|(b) 
shows that the strain to curvature ratio on the mid-line itself surpasses the classical buckling threshold value 
by nearly 20%. However, as Figure |l^(c) shows, this ratio drops away quickly in the direction transverse 
to the ridge-line. Thus the line plotted in Figure |l^(b) is a very localized maximum profile. The buckling 
envelope has a width on the order of the ridge width, and the strain to curvature ratio should be super- 
critical over this width before buckling occurs. The minimum cylinder buckling wavelength shrinks as A^^'^ 
as the sheet gets thinner, while the ridge width R shrinks as A^/^, so the cylinder buckling mode should be 
increasingly dependent only on local curvature and strain fields as A — > 0. Therefore the short wavelength 
localized mode should continue to be the preferred buckling mode for thinner sheets. Numerically the time 
required to compute eigenmodes grew very quickly as we decreased the thickness aspect ratio A, so we were 
not able to track the evolution of the short wavelength modes for significantly thinner sheets. 

4.4 Other geometries 
4.4.1 Different ridge angles 

Most of our in depth buckling mode data analysis has been performed on ridges with dihedral angle ^. 
However, the equation we derived for the critical stress (Eq. |5^) only depends on two parameters of the ridge 
shape - the transverse curvature and the material thickness. In order to show that this relation holds for 
more general geometries, we again consider ridges with different dihedral angles. As in Section 3.2.1, our 
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Figure 17: Curvature and strain profiles on ridge 
Plot (a) shows the curvature in units of ("+" symbol) and strain magnitude ("x" symbol) as a function 
of position along the ridge-line for a ridge at its buckling threshold. The ridge has a thickness aspect ratio 
of A = 5.5 X 10~*. Plot (b) shows the ratio of curvature to strain magnitude along the ridge- line for the 
same ridge at rest ("□" symbol) and at the buckling threshold ("o" symbol). The location and extent of 
the localized vibrational mode from Figure |l^(b) is highlighted. Plot (c) shows the ratio of curvature to 
strain magnitude across the middle of the ridge for the same ridge at rest ("□" symbol) and at the buckling 
threshold ("o" symbol). The horizontal lines in (b) and (c) are the classical buckling value of the strain to 
curvature ratio predicted by Eq Bq. 
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Figure 18: Curvature and strain profiles 
Plot (a) shows the ratio of strain magnitude to curvature along the ridge-line at the buckling threshold for 
five different dihedral angles and three different thicknesses. The ridges had dihedral angles ranging from ^ 
to Regardless of ridge angle, the values of the strain/curvature ratio tend to be grouped for ridges with 
the same thickness aspect ratio. In the top grouping, A — 1.25 x 10"'', in the middle A = 8 x 10~^, and for 
the bottom grouping A = 5 x 10~^. Plot (b) shows the buckling threshold values of the strain magnitude 
along the ridge-line for ridges with A — 1.25 x 10~^ and dihedral angles ranging from ^ (top line) to 22. 
(bottom line). The horizontal lines in (b) are the classical buckling values of the strain to curvature ratio 
predicted by Eq|58|for each thickness. 



29 



Figure 19: Shape of simulational grids 
The dashed hnes indi cate t he location of the simulated ridge-line. The grid in (b) was used for the simulations 
described in Section 4.4.2. The grid in (a) was used for all other ridge simulations. 
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Figure 20: Strain to curvature ratios for ridges with longer flanks 
The strain to curvature ratios on the ridge-line at the buckling threshold are shown for our typical ridges 
("o" symbols) and for ridges with longer flanks ("□" symbols). Thickness aspect ratios A shown here range 
from 1.25 x 10"'^ (top curves) to 5.6 x 10"^ (bottom curves). The horizontal lines are the classical buckling 
values of the strain to curvature ratio predicted by Eq Bq for each thickness. 



simulations used sheets with the same size and length to width ratio, but the location and orientation of the 
reflective planes for the edge boundary conditions were changed (the connection between reflective planes 
and dihedral angle is illustrated in Figure ||) . The energy of ridges with dihedral angle Tr — 2a scales as a^^^ , 
so adjusting this angle changes the curvature across the ridge signiflcantly. Still, as Figure |l^ shows for 
ridges with dihedral angles ranging from f to ^ that the ratio of strain to curvature along the ridge at the 
buckling threshold was the same for ridges with the same aspect ratio A. This is true despite the variation 
in the buckling strain by a factor of 2 between the largest and smallest angled ridges. This further afflrms 
the cylinder buckling hypothesis. 

This observation may also explain the insensitivity of the observed fractional change in energy between 
resting ridges and those at the buckling threshold. Another consequence of Lobkovsky's treatment in is 
that the longitudinal strain and transverse curvature on the ridge-line scale with the same power of a. Since 
our ridges always buckle at the same value of the ratio of strain to curvature, the identical scaling of these 
quantities with a implies an identical fractional change in their values (and therefore the total energy of 
the ridge) between the resting and buckling threshold states. We expect this behavior to extend to dihedral 
angle up to tt, though we haven't demonstrated it. 

4.4.2 Longer ridge flanks 

So far, we have simulated ridges on a grid whose edges formed a perfect square, as shown in Figure p^a). 
The ridge-line extended between two corners of this grid, so its simulated flanks were right triangles. The 
non- vertex corners of the grid were a distance X/2 from the center of the ridge-line, where X is the ridge 
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Figure 21: Buckling pattern on ridge with long flanks 

This ridge had a thickness aspect ratios A of 1.25 x 10~'^. Lighting and shading were chosen to emphasize 
physical features. 

length. As a variation on this ridge geometry, we also simulated ridges with flanks which were twice as 
long as those for the typical simulations, using the grid shown in Figure [lS|(b). On this grid the non- vertex 
corners were a distance X from the center of the ridge-line. The boundary conditions were again reflective 
planes which were oriented to give the ridge a ^ dihedral angle. 

We simulated ridges with thickness aspect ratio A ranging from 1.25 x 10"^ to 1.25 x lO"'^. For these 
ridges we once again found that the buckling threshold energy was approximately 20% greater than the 
resting ridge energy. As Figure ^ shows, we also found that the ratio of strain to curvature along the 
ridge- line at the buckling threshold was nearly equal to the threshold values for our typical ridges. 

Interestingly, the peaks in the strain-curvature ratio along the ridge-line were consistently closer to the 
vertices for this geometry than they were for all the other geometries we studied (this is visible in the profiles 
shown in Figure ^0|). This suggests that the locations of these peaks are determined by boundary conditions. 
The appearance of the peaks is probably an effect of the ridge pulling on its mirror image. Predictably, the 
change in location of this peak also causes the ridge to buckle closer to its vertices, as shown in Figure |2^. 
This is strong evidence that buckling occurs near the first localized patch on the ridge-line where the strain- 
curvature ratio is super-critical. When this pulling is absent, we might expect the buckling region to move 
towards the center. Such a case occurs in Lobkovsky et al.'s minimal ridge , where there is no applied stress 
at the boundaries, and buckling occurs at the center of the ridge. The change in buckling location is also 
consistent with the observed locations of new ridges on a buckled tetrahedron in j^. Since the tetrahedron 
ridges have shorter flanks than the cube we use, the new ridges may be expected to form closer to the center 
of the ridge on the tetrahedron. 

4.4.3 Shell buckling 

We have derived the behavior of the ridge buckling mode near the buckling transition as a function of the 
total longitudinal stress on the ridge-line, with as little reference as possible to where this stress comes from. 
At this point, we wish to address the role of the resting ridge stress in the buckling transition. As we 
showed in Section |^, even in its "resting" state the ridge has significant longitudinal stress - the ridge stores 
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Figure 22: Shell vs. ridge buckling energies 
This graph shows the diflerence between resting and buckhng energies for sheUs ("+" symbols) and ridges 
("x" symbols). The numerical scaling fits had exponents consistent with —1/3. 



a nonzero fraction its total elastic energy in the stress field. The amount of work required to buckle the 
ridge is therefore only that required to increase the longitudinal stress from the resting value to the buckling 
threshold value. It seems an odd occurrence that the energy required to make the ridge should also do part 
of the work required to break it. It is even more intriguing since the scaling of ridge stress and critical stress 
are identical - if this were not so then any ridge above or below (depending on the relative scaling) a critical 
length would buckle spontaneously.^ As it is, there are most likely some geometrical constraints placed on 
allowed ridge configurations purely by virtue of the fact that the ridge stress is naturally on the same scale 
as the buckling stress. 

To find out whether or not ridge stress weakens the ridge, we numerically studied how the buckling 
threshold changed when the resting ridge stress was removed. To do this, we first found the minimum 
energy configuration of a resting ridge for a given thickness aspect ratio A. We then redefined all the lengths 
and curvatures in the sheet such that the resting ridge configuration had zero stain and curvature, and thus 
zero resting energy. If we denote the strains and curvatures of the resting ridge as j°j and C°j respectively, 
we can write distortions away from this state as 

llj = Irj - l°j (72) 

C'ij = Cij — Cj°-, (73) 

where jij and Cij are computed as before. These primed quantities were substituted into the strain and 
curvature energy equations, Eqs. |^ and 0, to make the energy for the resting ridge configuration identically 
zero. We refer to sheets which have intrinsic curvature and non-flat metrics as a shells. 

As with the ridges before, we buckled the shells by imposing a gradually increasing hard wall potential at 
the vertices. We found that for any given thickness, it takes more work to buckle the shell than it does the 
corresponding ridge, though each started with exactly the same geometry (see Figure ^2|). Though the two 
systems evolve differently from their initial states under the applied load. Figure |2^(a) confirms that each 
buckles at nearly the same ratio of total stress to total (intrinsic plus extrinsic) curvature. This is just the 



result that would be predicted from the cylinder theory in Sec. 4.1, since that buckling mode depends only 
on the total value of the radius of curvature i? as a geometric quantity, without reference to its energetic 
cost. Figure |2^ shows that the strain to curvature ratios along the ridge-line are nearly the same for ridges 
and shells with identical thickness aspect ratios A. The figure also shows that the maximum of this ratios is 
near the same place on both ridges and shells, though it is more sharply peaked on shells. 

Also, Figure |2^(b) shows that although the shell cannot sustain the total amount of longitudinal stress 
that the ridge holds at the buckling threshold, it can sustain more additional stress, starting from the resting 

^Spontaneous buckling is commonly observed numerically due to lattice effects in very thin sheets. Any departure from the 
scaling laws due to lattice size in simulations or irregularities in real materials could be sufficient to push the required resting 
ridge stress over the buckling threshold. 
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Figure 23: Shell vs. ridge buckling configurations 
In each graph, "*" symbols denote values for ridges at rest, "x" symbols denote values for ridges at their 
buckling threshold, and "+" symbols denote values for shells at their buckling threshold. The numerical 
scaling fits in (a) had exponents consistent with 1 while those in (b) had exponents consistent with 2/3. 
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Figure 24: Strain to curvature ratios for shells 
The strain to curvature ratios on the ridge- line at the buckling threshold are shown for ridges ( "o" symbols) 
and shells ("□" symbols) with thickness ratios A ranging from 1.25 x 10~^ (top curves) to 2.5 x 10~* (bottom 
curves). The horizontal lines are the classical buckling values of the strain to curvature ratio predicted by 
Eq for each thickness. 
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ridge geometry, than the ridge before it buckles. The reason that the shell cannot sustain the same total 
stress is that it grows much flatter with applied load than does the ridge, so it reaches the critical stress to 
curvature ratio at a lower value of both these quantities. The difference in the evolution of the curvature 
with applied loading is due just to the fact that the cross-ridge curvature is determined by an energy balance, 
and the energetic term are much different for these two systems. 



4.5 Universality of buckling energy 

It is striking that for all the ridge geometries we studied, the buckling threshold energy was approximately 
20% greater than the resting energy. Our present ridge theory is not sufficient to fully explain this ratio, but 
we take this opportunity to speculate about how universal it may be. 

In Lobkovsky et al. stated that the bending and stretching energies on ridges should obey a virial 
relation, with total bending energy five times greater than total stretching. Assuming that both energies 
are only significant on the ridge-line, then the virial relation also extends to typical bending and stretching 
energy densities. Taking the largest terms from Eqs. ^ and 0, this gives: 



12(l-jy2) VV 

which reduces to 



V60 



From Eq ^ the classical value of the breaking strain for = 1/3 is 

7,, =« 0.61 hCyy. (76) 

In our simulations we observed that the bending energy doesn't change significantly as the ridge-line is 
compressed by the applied force. If the ridge curvature stayed the same and the ridge strain increased from 
its resting value up to the classic value, then the ratio of resting energy to buckling threshold energy would 
be (5 + 1):(5 -t- (0.61/0. 13)^. The breaking energy would be approximately 4.6 times the resting energy. 
However, Figure ^ shows that the mid-ridge curvature decreases as force is applied. Between resting and 
buckling, the curvature decreases by nearly a factor of two. The effect of this flattening of the ridge-line is 
to decrease the breaking strain by a factor of two and the breaking stain energy by a factor of four. The 
corresponding breaking energy is only 1.77 times greater than the resting energy. 

The discrepancy between the predicted factor 1.77 and the observed factor of 1.2 is understandable, given 
the simplicity of our approximations. In reality the distribution of strain and curvature on the ridge-line 
are not identical. In real ridges we may expect the local ratios to vary by factors of order unity depending 
on boundary conditions. However, we presently cannot explain why the local curvature at the center of the 
ridge- line drops by a factor of two while the total bending energy remains constant. This factor of two seems 
to be universal throughout our simulations, and we predict that the curvature won't drop by significantly 
larger fractions for different ridge boundary conditions. Still, a detailed understanding of this factor remains 
an open question. 



4.6 Discontinuity at buckling 

To complete our study of the buckling transition, we comment briefly on the post-buckled state and its 
rapid growth from the unbuckled state. As mentioned above, the equilibrium configuration immediately 
after buckling contains at least one large additional ridge. The additional ridges appear suddenly - when we 
first see them they are already as long as the unbroken ridge was wide. A significant change in the elastic 
energy accompanies the transition, as shown in Figure |2^. It is also notable that the post-buckled state bears 
little resemblance to the short wavelength buckling mode which we credit with causing the transition. Since 
the von Karman equations are highly non-linear, the growth of the buckled state quickly passes beyond the 
regime where it is well modeled by our linear stability analysis. We therefore presume that when the the 
non-linear terms start to become important, they favor further growth of the buckled state. The net result 
is an energetic avalanche into a completely different state. 
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Figure 25: Post-buckled energy 
This plot shows the total rescaled elastic energy of several simulated ridges as a function of rescaled inward 
vertex displacement just before and after the ridges buckled. The ridges had aspect ratios A ranging from 
5 X 10"* to 5 X 10~^. According to the scaling analysis in Section ^ the rescalings of (A/Ai)"^''"^ for energy 
and (A/Ai)^'^^ for vertex displacement would collapse all the lines on to one for perfect ridge scaling. The 
observed discrepancies between the pre-buckling part of the lines is small on the scale of the entire ridge 
evolution. In these graphs Ai = 5 x 10~*. 




Figure 26: Effect of buckling mode on R 
The solid line represents an unbuckled local patch of surface with dominant radius of curvature 7?. The 
dashed line shows the post-buckled surface with buckling mode amplitude C and wavelength (n- From 
Eq. the curvature at point A is approximately — -S- while the curvature at B is approximately + ^ . 
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As the buckling mode grows, it will begin to significantly perturb the pre-existing stress and curvature 
fields a and 1 jR. Figure ^ illustrates how the growth of the buckling mode perturbs the large transverse 
curvature. From Eq. the maxima for which the displacement w grows radially inward will decrease 
the local transverse radius of curvature by —CjC^, where C is the buckling mode amplitude and is its 
transverse wavelength. The maxima which grow radially outward will increase the local curvature by the 
same amount. We can also calculate the additional longitudinal stress due to the buckling mode itself from 
Eqs. H and |^ combined with Eq. 



The coefficient is always positive. Our frame is defined with positive normal displacements pointing 
inwards (downward and into the page in Figure |27| ) , so Eq. implies that there is additional compression 
at the points of maximum inward deflection of the buckling pattern, and matching extension at the points of 
maximum outward deflection (negative stress results from compression). Thus, as the buckling mode grows, 
the local ratio of strain to curvature becomes 



By the stability condition given in Eq. the inward growing maxima become more unstable to further 
growth while the outward ones become less unstable. 

From these simple arguments it is clear that the out-of-plane force balance changes dramatically as the 
buckling mode grows and non-linear terms become significant. Also, this reasoning indicates that the inward 
growing maxima are more favorable than outward growing ones. On a cylinder, the constraints of periodicity 
require an equal azimuthal number of inward and outward maxima, so force balance is again achieved with 
the same number of maxima as the initial unstable mode. This constraint does not hold on the ridge, so 
inward maxima are free to grow into the region of the ridge flanks. The net result could be that one inward 
peak grows until it subsumes all the other maxima and becomes the one prominent feature of the buckled 
state - a single large transverse ridge. 

The only reason for the growth of the new ridge to cease is that the potential energy driving this motion 
is exhausted. Figure |^ shows how the elastic energy redistributes itself upon the buckling of a ridge with 
aspect ratio 10~^. Predictably, the largest local decrease in elastic energy is the loss of stress energy density 
along the ridge. In our buckling scheme it is the stored longitudinal stress which drives the growth of the 
buckling mode - we observed that the inward buckling peak which becomes the observed single new ridge 
grows until the ridge-line stress is nearly gone. Some of the energy is stored in the new ridge and new 
vertices, but less than was stored in the threshold state. 

Besides showing the jump in energy at buckling. Figure ^ also shows the randomness in the observed 
energies of the postbuckled state. Our simulations were optimized for the pre-buckled ridges. Their accuracy 
in modeling the post-buckled state is qualitative at best. There is no indication from this graph that the 
post-buckled state has the same energy scaling as the ridge. This is not surprising, since the "resting" 
configuration of our new post-buckled state is not defined, so we cannot be sure that the states we see 
immediately after buckling are at "equivalent" values of the tip displacement, in the sense developed in 
Section ^ for a similarity solution at different material thicknesses. 

Some excellent experimental work on the preferred crease size for a circular cylindrical cross section under 
axial compression is presented in |l^ . In this work the authors found that the actual saturation length 
of the new crease was determined by a balance between the energy of the crease and that of the additional 
singularities at its ends. In our geometry, the scaling of the new ridge energy changes once it has grown 
beyond the width of the original boundary layer. Therefore, the energy balance determined for the constant 
curvature cross section may not be applicable. We leave this topic for future research. 
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Figure 27: Redistribution of elastic energy at buckling 
Each plot shows the postbuckling configuration first shown in Figure |l^, with one large new ridge crossing 
the original ridge. The original ridge had aspect ratio A — 10~^. The plots are shaded according to change in 
local elastic energy density between the buckling threshold values and the postbuckled value for the pictured 
configuration - no change is gray, increases are white, and decreases are black. Image (a) shows the change 
in total elastic energy, image (b) shows the change in bending energy, and image (c) shows the change in 
stretching energy. 
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5 Discussion 



We have explored the behavior of a stretching ridge under the apphcation of an external force potential. 
For the sake of clarity we have focused our simulations on a particular representative ridge geometry, but 
the response of this ridge has been shown to obey very general principles. In this section we re-cap our 
discoveries, putting them into the broader context of the enhanced strength these spontaneous structures 
add to thin sheets. We also discuss the range of applicability for our approach to other ridge geometries, 
and to the behavior of collections of ridges in a crumpled sheet. Finally, we suggest engineering applications 
of the understanding we have gained concerning stretching ridges. 

The resistance of materials to typical forms of distortion and damage is a very well established field, with 
a history that dates back to the 19th century. However, a crumpled sheet derives its strength not just from 
its material properties, but also from the spontaneous ridge network it contains. This spontaneous network 
confers strength in a way that clearly arises from the co-operative interaction between curvature and strain. 
The novel aspect of this interaction has already been shown by the newly identified scaling of the energy of 
these structures with overall size of the system However, the strength against collapse resulting from 
these structures has up to now been poorly understood. In a highly crumpled sheet, resistance to further 
deformation results almost entirely from the work required to deform and break the ridges which span the 
volume occupied by the sheet. The strength of ridges in turn results from their shape, and their effective 
elastic modulus is not related to the modulus of the component material in any simple way. 

In Sections ^ and || we established a scaling relation for the response of a ridge to forces applied at its 
endpoints. This is the type of forcing against which ridges are strongest (have the highest effective modulus). 
Presumably, when a force is applied in an arbitrary direction to a moderately crumpled sheet, ridges which 
are oriented at broad angles to the applied force will yield very quickly to it, and resistance to the force will 
come from ridges which happen to be aligned parallel to the forcing. Thus the ridge response to this particular 
forcing determines the effective elastic modulus of crumpled sheets. We showed that, given a knowledge of 
how one ridge of any size will respond to the force at its ends, we can rescale the force to displacement 
relation to all other ridge lengths by multiplying it by a simple power of the thickness aspect ratio A, namely 
A"'^. The force to displacement relation for an individual ridge can be obtained through simulations or simple 
estimates. Together with a model of the distribution of ridge sizes in a typical crumpled sheet, our scaling 
relation for the ridge strength gives a complete model for the effective elastic modulus of the entire crumpled 
sheet, as well as the change in the sheet's strength as it is further crumpled and the typical ridge sizes change. 

We derived the scaling of the ridge force response by first assuming that ridge scaling was still valid for 
forced ridges, then calculating the required rescaling of the perturbing force. Our approach is limited by the 
requirement that both the location and magnitude of applied forcing must be rescaled for a similarity solution. 
Also, there is no systematic way to determine if the scaling assumption will hold for every forcing. Still, our 
approach is comparable to other treatments in its effectiveness, since much of the physics of crumpling relies 
on intuition for each special case. 

Happily, our approach was shown to be well suited for point forces applied to the vertices. For this 
case, the location of the forcing is fixed under rescaling, since the vertices by definition do not move when 
the ridge gets thinner. Also, since the applied forcing works almost entirely to compress the ridge-line, its 
coupling to the longitudinal ridge stress will be very strong. Thus the predicted scaling of the force response 
is unambiguous. 

Forcing applied to other points on the boundary of the sheet will most likely not scale as cleanly as forcing 
applied to the vertex. For other locations and angles, the applied force may result in large stress transverse 
to the ridge - it would therefore strongly perturb both transverse and longitudinal stresses. However, these 
stresses have different scalings on the resting ridge, so an equal perturbation of each would most likely ruin 
the ridge scaling. 

The other important perturbation, which we did not simulate, is forcing applied normal to the surface. 
Scaling of this force response, as derived in Section ^, should be fairly robust since the term P is perturbing 
a quantity which is zero for resting ridges. Therefore there is no pre-established scaling to destroy. 

Our other important result, established in Section ^, was to link the buckling transition for ridges to 
the buckling of thin cylinders. This result is supported by a great deal of analysis and is very general. 
Prior speculation held that ridges may derive anomalously large breaking strength from their pre-existing 
longitudinal strain. We have shown that, in terms of strength, the ridge acts essentially as a cylinder 
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whose radius scales with thickness. Whereas the work in Sections || and ^ allowed us to understand the 
strength of crumpled sheets against small deformations which did not change the structure of the crumpling 
network, knowledge of the buckling strength lets us model the evolution of the strength and energy of a 
sheet throughout the crumpling process, from the flat to the highly crumpled state. 

In Section |4[ we show that the stability of the ridge against buckling is determined completely by the 
local ratio of the transverse strain to the longitudinal curvature on the ridge-line. Since the transverse 
curvature on the ridge scales with its length, we can immediately determine the buckling stress of any ridge 
as a function only of its length and thickness. We established that the allowed buckling wavelengths Ccj are 
between 2ttX\^/^ < C,ci < X, where X is the ridge length and A is the thickness aspect ratio. Buckling 
can take place when the strain to curvature ratio is supercritical over a region larger than the minimum 
wavelength. We showed that ridges buckle near the point at which this ratio has a localized maximum on 
the ridge-line. We established that the location of this strain to curvature maximum does not depend on 
the dihedral angle of the ridge, but it does depend on the angle of the ridge-line relative to its neighboring 
ridges. The simplicity of the buckling criterion established here, as well as the clear connection between this 
ratio and the buckling behavior of ridges, is a great improvement over the previous understanding of the 
breaking strength of these structures. Further development of the relation between a ridge and its neighbors 
may lead to general laws regarding preferred distributions of angles separating ridges in crumpled sheets - 
this plus the length and energy distributions discussed below could lead to an accurate statistical mechanics 
for crumpled sheets 

It is our hope that the knowledge gained in this study can be helpful in the development of a statistical 
mechanics for ridge distributions in crumpled sheets. We have demonstrated for a range of ridge angles under 
a typical form of forcing that the energy at the buckling threshold is a fixed multiple of the resting ridge 
energy. For our measurements this multiple was approximately 1.2 - we speculate based on the arguments 



in Section 4.5 that this multiple will not be greater than 2 in general. Combined with previous work which 
uncovered the length and angle scaling of the ridge energy, this limitation on the possible loading of a 
ridge should place some limitation on the relative sizes of adjacent ridges in a typical crumpled sheet. This 
proposed limitation follows from the assumption that in the interior of the sheet, much of the force on 
vertices is carried through adjacent ridges. Knowledge of the strength of ridges will definitely lend insight 
to the evolution of successive crumpled states as a sheet is compressed. 

Finally, in terms of applicability to real-world problems, the understanding of ridge buckling developed 
here has practical import for the possible use of single ridges as structural elements. We observed that the 
weakest point on the ridge is near the point of largest stress to curvature ratio. Thus ridges which are used as 
support elements could be reinforced selectively at areas determined to be weak points through this analysis. 
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